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Abstract. In this article, wc address both recent advances and open questions 
in some mathematical and computational issues in geophysical fluid dynamics 
(GFD) and climate dynamics. The main focus is on 1) the primitive equa- 
tions (PEs) models and their related mathematical and computational issues, 
2) climate variability, predictability and successive bifurcation, and 3) a new 
dynamical systems theory and its applications to GFD and climate dynamics. 
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1. Introduction 

The atmosphere and ocean around the earth are rotating geophysical fluids, 
which are also two important components of the climate system. The phenomena of 
the atmosphere and ocean are extremely rich in their organization and complexity, 
and a lot of them cannot be produced by laboratory experiments. The atmosphere 
or the ocean or the couple atmosphere and ocean can be viewed as an initial and 
boundary value problem (Bjerknes [5], Rossby |125j . Phillips [120j ). or an infinite 
dimensional dynamical system. These phenomena involve a broad range of temporal 
and spatial scales (Charney [H]). For example, according to J. von Neumann |147j . 
the motion of the atmosphere can be divided into three categories depending on 
the time scale of the prediction. They are motions corresponding respectively to 
the short time, medium range and long term behavior of the atmosphere. The 
understanding of these complicated and scientific issues necessitate a joint effort 
of scientists in many fields. Also, as John von Neumann [147| pointed out, this 
difficult problem involves a combination of modeling, mathematical theory and 
scientific computing. 

In this article, we shall address mathematical and numerical issues in geophysical 
fiuid dynamics and climate dynamics. The main topics include 

(1) issues on the modeling, mathematical analysis and numerical analysis of 
the primitive equation (PEs), 

(2) climate variability, predictability and successive bifurcation, 

(3) a new dynamical systems theory and its applications to geophysical fluid 
dynamics. 

As we know, the atmosphere is a compressible fluid and the seawater is a slightly 
compressible fluid. The governing equations for either the atmosphere, or the ocean, 
or the coupled atmosphere-ocean models are the general equations of hydrodynamic 
equations together with other conservation laws for such quantities as the energy, 
humidity and salinity, and with proper boundary and interface conditions. Most 
general circulation models (GCMs) are based on the PEs, which are derived using 
the hydrostatic assumption in the vertical direction. This assumption is due to the 
smallness of the aspect ratio (between the vertical and horizontal length scales). We 
shall present a brief survey on recent theoretical and computational developments 
and future studies of the PEs. 

One of the primary goals in climate dynamics is to document, through careful 
theoretical and numerical studies, the presence of climate low frequency variability, 
to verify the robustness of this variability's characteristics to changes in model pa- 
rameters, and to help explain its physical mechanisms. The thorough understanding 
of this variability is a challenging problem with important practical implications 
for geophysical efforts to quantify predictability, analyze error growth in dynami- 
cal models, and develop efficient forecast methods. As examples, we discuss a few 
sources of variability, including wind-driven (horizontal) and thermohaline (vertical) 
circulations. El Niho-Southern Oscillation (ENSO). and Intraseasonal oscillations 
(ISO). 

The study of the above geophysical problems involves on the one hand applica- 
tions of the existing mathematical and computational theories to the understanding 
of the underlying physical problems, and on the other hand the development of new 
mathematical theories. 
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We shall present briefly a dynamic bifurcation and stability theory and its appli- 
cations to GFD. This theory, developed recently by Ma and Wang |100j . is for both 
finite and infinite dimensional dynamical systems, and is centered at a new notion 
of bifurcation, called attractor bifurcation. The theory is briefly described by a 
simple system of two ordinary differential equations, and by the classical Raylcigh- 
Benard convection. Applications to the stratified Boussinesq equations model and 
the doubly-diffusive models are also addressed. 

We would like to mention that there are many important issues not covered in 
this article, including, for example, the ocean and atmosphere data assimilation 
and prediction problems, and the stochastic-dynamics studies; see, among many 
others, [SlllSHlEOllMlESllMlESlEIllMllSTlE^ \M\ and the references 

therein. 

The article is organized as follows. In Section 2, some basic GFD models are 
introduced, with some mathematical and computational issues given in Section 3. 
Section 4 is on predictability, and Section 5 deals with issues on climate variability. 
Section 6 presents the new dynamical systems theory based on attractor bifurcation 
and its application to Rayleigh-Bcnard convection and to GFD models. 

2. Modeling 

2.1. The primitive equations (PEs) of the atmosphere. Physical laws gov- 
erning the motion and states of the atmosphere and ocean can be described by 
the general equations of hydrodynamics and thermodynamics. Using a non-inertial 
coordinate system rotating with the earth, these equations can be written as follows: 

dV ^ ^ ^ ^ 1 
— + y • V3F + 2nxV-g+ -gradap = Dm, 
ot p 

^+diy,ipV)^0, 

^^■^^ Cv^+c^V-VsT+^V-iV^Q + DH, 

ot p 



ot p 
p = RpT. 

Here the first equation is the momentum equation, the second is the continuity 
equation, the third is first law of thermodynamics, the fourth is diffusion equation 
for the humidity, and the last is the equation of state for ideal gas. The unknown 
functions are the three dimensional velocity field V, the density function p, the 
pressure function p, the temperature function T, and the specific humidity function 
q. Moreover, in the above equations, n stands for the angular velocity of the 
earth, g the gravity, R the gas constant, Cy the specific heat at constant volume. 
Dm the viscosity terms, Dh the temperature diffusion, Q the heat flux per unit 
density at the unit time interval, which includes molecule or turbulent, radiative 
and evaporative heating, and S the differences of the rates of the evaporation and 
condensation. 

These equations are normally far too complicated; simplifications from both the 
physical and mathematical points of view are necessary. There are essentially two 
characteristics of both the atmosphere and ocean, which are used in simplifying the 
equations. The first one is that for large scale geophysical flows, the ratio between 
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the vertical and horizontal scales is very small; this leads to the primitive equations 
(PEs) of both the atmosphere and the ocean, which are the basic equations for 
these two fluids. More precisely, the PEs are obtained from the general equations 
of hydrodynamics and thermodynamics of the compressible atmosphere, by ap- 
proximating the momentum equation in the vertical direction with the hydrostatic 
equation: 

(2.2) = -pg. 

OZ 

This hydrostatic equation is based on the ratio between the vertical and horizontal 
scale being small. Here p is the density, g the gravitational constant, and z = r — a 
height above the sea level, r the radial distance, and a the mean radius of the 
earth. Equation (|2.2p expresses the fact that p is a decreasing function along the 
vertical so that one can use p instead of z as the vertical variable. Motivated by 
this hydrostatic approximation, we can introduce a generalized vertical coordinate 
system s-system given by 

(2.3) s = s{e,^,z,t), 



where s is a strict monotonic function of z. Then the basic equations of the large- 
scale atmospheric motion in the s-system are 



(2.4) 



d'>^ r-, .dv ^, 1 „ ^ 

— + u • VsU + s— + fk XV + -VsZ = Dm, 

Ot OS p 

dpds 

-^^+P9 ^ 0, 
OS OZ 



d f dp 
d's \ dt 
dT 



V . ( v—"] + — ( li- 
ds J ds \ ds 



xirra. -dT I dp dp 
dt ds p \ dt ds 

9q „ .dq S ^ 
-j+w V.g + STT = -+Dg. 
dt ds p 



= Q + Dh, 



Some common s-systems in meteorology are respectively the p-system (the pres- 
sure coordinate), the a-system (the transformed pressure coordinate), the 6'-system 
(the isentropic coordinate), and the (^-system (the topographic coordinate or trans- 
formed height coordinate). The above PEs appear in the literature in e.g. the 
books of A. E. GiU [42], G. Haltiner and R. WiUiams [43], J. R. Holten [44], J. 
Pedlosky [TIH], J. P. Peixoto and A. H. Oort [Hll, W. M. Washington and C. L. 
Parkinson |152| . Q. C. Zeng [159j . We remark here that sometimes the pressure 
coordinate is denoted by ij, and terrain-following by a. 

For simplicity, here we discuss only the case with the coordinate transformation 
from {9^ip,z) to {9,ip,p). The basic equations of the atmosphere are then the 
Primitive Equations (PEs) of the atmosphere in the p-coordinate system. As they 
appear in classical meteorology books (see e.g. [159] and Salby [128j ). the PEs are 
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given by 



dv dv 

— + v ■ ^v + uj— + 2ilcos6'fc X u + V$ = Dm, 
at op 



0, 



(2.5) div V 



9$ RT 
dp p 

dp 

— + v -VT + LU I — = \ h Dh, 

Ot \ p dp J Cp Cp 



^+vWq + LO^ = E-C + Dg, 

where Dm is the dissipation terms for momentum and Dh and Dq are diffusion 
terms for heat and moisture, respectively, E and C are the rates of evaporation and 
condensation due to cloud processes, Cp the heat capacity, and Qrad and Qcon the 
net radiative heating and the heating due to condensation processes, respectively. 
We use the pressure coordinate system {d,Lp,p), where 9 {0 < 9 < tt) and f {0 < 
if < 2-k) are the colatitude and longitude variables, and p the pressure of the 
air. The nondynamical processes Qrad, Qcon, E and C are called model physics. 
Furthermore, the unknown functions are the horizontal velocity v, the vertical 
velocity lo = dp/dt, the geopotential $, the temperature T, and specific humidity 
q. The operators div and V are the two dimensional operators on the sphere. 

Of course, this set of equations is supplemented with a set of physically sound 
boundary conditions such as p.3p . depending on the specific form of the forcing 
and dissipation. 

2.2. Ocean models. The sea water is almost an incompressible fluid, leading to 
the Boussinesq approximation, i.e., a variable density is only recognized in the 
buoyancy term and the equation of state. The resulting equations are called the 
Boussinesq equations given as follows: 

dv dv 1 d'^v 

— + VyV + w- — I gradp + jk x v — fiAv - v^r-^ = 0, 

dt dz po dz'^ 



dw dw 1 dp p d^w 

— + VyW + w— H — H g - pAw ~ v—^ = 0, 

dt dz Po dz Po dz'' 

(2.6) = 

dT dT , d^T 

— + V,T + w— - pT^T - = 0, 
dt dz dz'^ 

dS ^ ^ dS ^ ^ d^S ^ 

— + VyS + w- /isAA - i^s^-Y = 0' 

dt dz dz'' 

p = po{l - I3t{T - To) + Ps{S - So)), 

where v is the horizontal velocity field, w the vertical velocity, and S the salinity. 
The sixth equation in (j2.6|l is an empirical equation for the density function based 
on the linear approximation. In general, density p is a nonlinear function of T, S, 
and p. With higher approximations, one will encounter additional mathematical 
difficulties although the nonlinear equation of state is essential for some elements 
of ocean circulation (e.g., cabbehng). 
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As in the atmospheric case, the hydrostatic assumption is usuaUy used, leading 
to the PEs for the large-scale ocean: 

dv dv 1 d'^v 

— + VvV + W— H \/p + jk X V ~ IMviXaV - Vv^r^ = 0, 

at oz po az'^ 

dp 

oz 

(2.7) " + 

dT dT ^ d^T 

— + V„T + W— HT^aT - vt-^ = 0, 

at oz oz^ 

dS ^ ^ dS ^ „ d^S ^ 

— + V„6 + W- /XsAaA - = 0, 

at oz oz^ 

p = po{l - Pt{T - To) + l3siS^ So)). 

Also, we note that if the hydrostatic assumption is made first, the Boussinesq 
approximation is not really necessary (see eq. (2.5) - divergcncc-frec! - or, e.g., de 
Szoeke and Samelson [TO] . 

2.3. Coupled atmosphere-ocean models. Coupled atmosphere and ocean mod- 
els consist of 1) models for the atmosphere component, 2) models for the ocean 
component, and 3) interface conditions. The interface conditions are used to cou- 
ple the atmosphere and ocean systems, and are usually derived based on on first 
principles; see Gill [12], Washington and Parkinson |152| . A mathematically well- 
posed coupled model with physically sound interface conditions and the PEs of the 
atmosphere and the ocean are given in Lions, Temam and Wang [55] • We refer 
interested readers to these references to further studies. 

As we know, the atmosphere and ocean components have quite different time 
scales, leading to complicated dynamics. For example, from the computational 
point of view, one needs to incorporate the two time scales; see e.g. }88) . 



3. Some Theoretical and Computational Issues for the PEs 

3.1. Dynamical systems perspective of the models. From the mathematical 
point of view, we can put the models addressed in the previous section in the 
perspective of infinite dimensional dynamical systems as follows: 

ipt + + R{ip) = F, 
(3.1) . 

<<5|t=0 = </'o, 

defined on an infinite dimensional phase space H. Here A : H H is an unbounded 
linear operator, R : H ^ H is a nonlinear operator, F is the forcing term, and ipo 
is the initial data. 

We remark here that the linear operator can usually be written as A = Ai + A2, 
where Ai stands for the irreversible diabatic linear processes of energy dissipa- 
tion, and A2 for the reversible adiabatic linear processes of energy conversation. 
The nonlinear term R{(p) represents the reversible adiabatic nonlinear processes of 
energy conversation. The properties of these operators refiect directly the essen- 
tial characteristics of two kinds of basic processes with entirely different physical 
meanings. 
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The above formulation is often achieved by (a) establishing a proper functional 
setting of the model, and (b) proving the existence and uniqueness of the solutions. 

Hereafter we demonstrate the procedure with the PEs. Due to some technical 
reasons, some minor and physically reasonable modifications of the PEs are made. 
In particular, we assume that the model physics Qrad, Qcon, E and C are given 
functions of location and time. We specify also the viscosity, diffusion terms as (see 
among others Lions et. al. [53] and Chou [H]): 

Dm = -Liv, 

Qrad Qcon „ r i /O 
1 h Uh = -L2I + Qt, 



E-C + Dg = -L3q + Qq, 

where /i^, Vi are horizontal and vertical viscosity and diffusion coefficients, A is the 
Laplace operator on the sphere, Qt and Qq are treated as given functions, and 
T{p) a given temperature profile, which can be considered as the climate average 
of T. The boundary conditions for the PEs are given by 

d 

T, q) = {js{vs - v),as{Ts - T), aq{qs - q)}, u; = at p = P, 

(3.3) 

■Q;^{v,T,q) ^ Q, uj ^ at p = po. 

The second and third equations (j2.5|l are diagnostic ones; integrating them in 
p-direction, we obtain 

fP 

div v{p )dp = 0, 

'Po 

(3.4) uj = W{v) = - / div v{p')dp', 

Jpo 
J p P 

Then the PEs arc equivalent to the fohowing functional formulation: 
Ou 

— + A{v)u + P{u) + Lu+ (V$„ 0, 0) = (0, Qt, Qg), 

(3.5) p 

div / vdp = 0, 

du 

where u = (v, T, q), A(v)u — VyU + W(v) — , Lu corresponds to the viscosity and 

dp 

diffusion terms, and Pu the lower order terms. The above new formulation was 
first introduced by Lions, Temam and Wang [83] . 

We solve then the PEs in some infinite dimensional phase spaces H and V . In 
particular, we use 



Hi = \v(^L^ I div I vdz = 
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as the phase space for the horizontal velocity v. Then we project in the phase 
space. Using this projection, the unknown function <i>s plays a role as a Lagrangian 
multiplier, which can be recovered by the following decomposition: 



Then the PEs are equivalent to an infinite dimensional dynamical system as (|3.ip . 

With the above formulation, for example, we encounter the following new non- 
local Stokes problem: 



supplemented with suitable boundary conditions. From the mathematical point of 
view, all techniques for the regularity of solutions are local. But our problem here 
is nonlocal; the regularity of the solutions for this problem can be obtained using 
Nirenberg's finite difference quotient method; see [83] . 

Other models such as the PEs of the ocean and the coupled atmosphere-ocean 
models can be viewed as infinite dimensional dynamical systems in the form of 
p.ip . We refer the interested readers to Lions, Temam and Wang [84] for the PEs 
of the ocean, and |85l [87] for the coupled atmosphere-ocean models. 

3.2. Well-posedness. One of the first mathematical questions is the existence, 
uniqueness and regularity of solutions of the models. The main results in this 
direction can be briefly summarized as follows, and we refer the interested readers 
to the related references given below for more details: 

For the PEs of the large-scale atmosphere, the existence of global (in time) weak 
solutions of the primitive equations for the atmosphere is obtained by Lions, Temam 
and Wang [53], where the new formulation described above is introduced. 

In fact, the key ingredient for most, if not all, existence results for the PEs 
of the atmosphere-only, the ocean-only, or the coupled atmosphere-ocean depend 
heavily on the formulation p.4p and (|3.5p . first introduced by Lions, Temam and 
Wang [55]. We note that without using this new formulation, one can also obtain 
the existence of global weak solutions by introducing proper function spaces with 
more regularity in the p-direction; see Wang [150] . However, the new formulation 
is important for viewing the PEs as an infinite dimensional system. 

The existence of global strong solutions is first obtained for small data and large 
time or for short time by [150|, 1141] . Other studies include the case for thin domain 
case HH] (see also Ifitimie and Raugel |5D] for discussions on the Navier-Stokes 
equations on thin domains), and for fast rotation [1]. 

Recently, for the primitive equations of the ocean with free top and bottom 
boundary conditions, an existence of long time strong solutions with general data is 
obtained by Cao and Titi [5] , using the new formulation and the fact that the surface 
pressure depends on the two spatial directions, and then similar result is obtained 
for the Dirichlet boundary conditions by Kukavica and Ziane |65| . Furthermore, Ju 
[54] studied the global attractor for the primitive equations. 

The corresponding results can also be obtained for the PEs of the ocean; see 
[84l I141|, [8] and the references therein for more details. Existence of global weak 
solutions were obtained for the coupled atmosphere-ocean model introduced in [87] . 



Av + = /, 



(3.6) 




Po 
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As mentioned earlier, the hydrostatic assumption and its related formulation 
given by (|3.5|) and ()3.6|) are crucial in many of the existence results for both strong 
and weak solutions. We would like to point out that Laprise [B5j suggests that 
hydrostatic-pressure coordinates could be used advantageously in nonhydrostatic 
atmospheric models based on the fully compressible equations. We believe that 
with the Laprise' formulation, one can extend some of the results discussed here to 
non-hydrostatic cases. 

Another situation where the new formulation (|3.5p and (j3.6p does not appear 
to be available is related to more complex low boundary conditions. Some partial 
results for weak solutions are obtained in |154| , and apparently many related issues 
are still open. 

3.3. Long-time dynamics and nonlinear adjustment process. Regarding the 
PEs as an infinite dimensional dynamical system, the existence and finite dimen- 
sionality of the global attractors of the PEs with vertical diffusion were explored in 
[83l EH [82] ■ The finite dimensionality of the global attractor of the PEs provides a 
mathematical foundation that the infinite dynamical system can be described by a 
finite dimensional dynamical system. 

As we know all general circulation models of either the atmosphere-only, or the 
ocean-only, or the coupled ocean-atmosphere systems are based on the PEs with 
more detailed model physics. In both these GCM models and theoretical climate 
studies, the PEs are often replaced by a set of truncated ordinary differential equa- 
tions (ODEs), whose asymptotic solution sets, called attractors, can be investigated 
in the more tractable setting of a finite-dimensional phase space, without seriously 
altering the essential dynamics. It is, however, not known mathematically whether 
this truncation is really reasonable, and moreover, how we can determine a priori 
which finite-dimensional truncations are sufficient to capture the essential features 
of the atmosphere or the oceans. With this objective, nonlinear adjustment process 
associated with the long-time dynamics of the models were carefully conducted in 
a series of papers by Li and Chou l7Tl[72l[23l[7l[75l[76l[7Zl[78];see also [mj [82] . 

An important consequence of the above mentioned results for the long-time dy- 
namics is that the nonlinear adjustment process of the climate is a forced, dissi- 
pative, nonlinear system to external forcing. This nonlinear adjustment process is 
different from the adjustments in the traditional dynamical meteorology, including 
the geostrophic adjustment, rotational adjustment, potential vorticity adjustment, 
static adjustment, etc. 

These traditional adjustments do not appear to be associated with any attract- 
ing properties that attractors usually process. It is indicated from the nonlinear 
adjustment process that as time increases, the information carried by the initial 
state will gradually be lost. In addition, there are three categories of time bound- 
ary layers (TBL) and the self-similar structure of the TBL in the adjustment and 
evolution processes of the forced, dissipative, nonlinear system. 

An important question is how to determine the structure of the attractors and 
the distribution of their attractive domains under the given conditions of known 
external forcing. Although the information of initial state will be decayed as the 
time increases, it does not mean that the information of initial state is not important 
to long-time dynamics. The quantity of those initial values which locate very close 
to the boundary between two different attractive domains are very important since 
their local asymptotic behavior will be quite different due to slight initial error. 
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Another open question is how to construct independent orthogonal basis of the 
attractor in practice based on the finite dimensionahty of the global attractor. Two 
empirical approaches used at present are the Principal Analysis (PC) or Empirical 
Orthogonal Function (EOF) and Singular Vector Decomposition (SVD) based on 
the time series of numerical solutions or observational filed data. They are available 
although they are empirical. Qiu and Chou [122| and Wang et al [149j made a 
valuable attempt to apply the long-time dynamics of the atmosphere mentioned 
above and this two empirical decomposition methods for orthogonal bases of the 
attractor to study 4-dimensional data assimilation. 

3.4. Multiscale asymptotics and simplified models. As practiced by the ear- 
lier workers in this field such as J. Charney and John von Neumann, and from 
the lessons learned by the failure of Richardson's pioneering work, one tries to be 
satisfied with simplified models approximating the actual motions to a greater or 
lesser degree instead of attempting to deal with the atmosphere in all its complex- 
ity. By starting with models incorporating only what are thought to be the most 
important of atmospheric influences, and by gradually bringing in others, one is 
able to proceed inductively and thereby to avoid the pitfalls inevitably encountered 
when a great many poorly understood factors are introduced all at once. 

One of the dominant features of both the atmosphere and the ocean is the influ- 
ence of the rotation of the earth. The emphasis on the importance of the rotation 
effects of the earth and their study should be traced back to the work of P. Laplace 
in the eighteenth century. The question of how a fluid adjusts in a uniformly ro- 
tating system was not completely discussed until the time of C. G. Rossby |126j . 
when Rossby considered the process of adjustment to the geostrophic equilibrium. 
This process is now referred to as the Rossby adjustment. Roughly speaking the 
Rossby adjustment process explains why the atmosphere and ocean are always close 
to geostrophic equilibrium, for if any force tries to upset such an equilibrium, the 
gravitational restoring force quickly restores a near geostrophic equilibrium. Later 
in the 1940's, under the famous /3-plane assumption, Charney [13] introduced the 
quasi-geostrophic (QG) equations for the large-scale (with horizontal scale compa- 
rable to 1000 km) mid-latitude atmosphere. Since then, there have been many 
studies from both the physical and numerical points of view. This model has been 
the main driving force of the much development of theoretical meteorology and 
oceanography. 

Mathematically speaking, the QG theory is based on asymptotics in terms of a 
small parameter, called Rossby number 7?o, a dimensionless number relating the 
ratio of inertial force to Coriolis force for a given flow of a rotating fluid. The key 
idea in the geostrophic asymptotics, leading to the QG equations, is to approximate 
the spherical midlatitude region by the tangent plane, called the /?-plane, at the 
center of the region, and to express the Coriolis parameter in terms of the Rossby 
number. 

Thanks in particular to the vision and effort of Professor J. L. Lions, there have 
been extensive mathematical studies. As this topic is very well-received and studied 
by applied mathematicians, we do not go into details, and the interested readers are 
refereed to (Lions, Temam & Wang [86l[89j, Bourgeois & Beale [6], Babin, Mahalov 
& Nicolaenko [1], Gallagher [29], Embid & Majda [25]) and the references therein. 
In addition, planetary geostrophic equations (PGEs) of ocean have also received a 
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lot of attention recently following the early work by Samelson, Temam and Wang 
[13H I130j . However, many issues are still open. 

As mentioned earlier, many geophysical processes have multiscale characteristics. 
One aspect of the studies requires a careful examination of the interactions of 
multiple temporal and spatial scales. A combination of rigorous mathematics and 
physical modeling together with scientific computing appear to be crucial for the 
understanding of these multiscale physical processes. ENSO and ISO are two such 
examples (see also Section 5). 



3.5. Some computational issues. On the one hand, we need to develop more 
efficient numerical methods for general circulation models (GCMs), including at- 
mospheric general circulation model (AGCM), oceanic general circulation model 
(OGCM) and coupled general circulation model (CGCM), climate system model 
and earth system model. On the other hand, numerical simulations are used to 
test the theoretical results obtained, as well as for preliminary exploration of the 
phenomena apparent in the governing PDEs, and to obtain guidance on the most 
interesting directions for theoretical studies. 

Here we simply address some computational issues without discussing physical 
processes, which are certainly important in developing GCMs. Two basic discretiza- 
tion schemes commonly used in GCMs are grid-point and spectral approaches. 
There are many crucial issues which are not fully resolved, including 1) spherical 
geometry and singularity near the polar regions, 2) irregular domains, 3) multiscale 
(spatial and temporal) problems involving both fast and slow processes, 4) vertical 
stratification, 5) sub-grid processes, 6) model bias, and 7) nonhydrostatic mod- 
els, etc. We note that although many existing models are formulated and solved in 
spherical coordinates, numerical and mathematical difficulties caused by the depen- 
dence of the meshes size on the latitude is still not fully resolved. Furthermore, high 
precision computation and very fine resolutions are two tendencies in developing 
numerical models. 

Recently a fast and efficient spectral method for the PEs of the atmosphere is 
introduced by Shen and Wang |133| , and further studies for this method applied to 
more practical GCMs are needed. Another paper uses heavily the surface pressure 
formulation of the PEs for the ocean is given in Samelson et al. [129j . Its incor- 
poration to OGCMs, and its simulation for studying specific oceanic phenomena 
appear to be necessary. 

A remarkable, but neglected, problem is on influences of round-off error on long 
time numerical integrations since round-off error can cause numerical uncertainty. 
Owing to the inherent relationship between the two uncertainties due to numerical 
method and finite precision of computer respectively, a computational uncertainty 
principle (CUP) could be definitely existed in numerical nonlinear systems (Li et 
al. |80l 181] ) , which implies a certain limitation to the computational capacity of 
numerical methods under the inherent property of finite machine precision. In prac- 
tice, how to define the optimal time step size and optimal horizontal and vertical 
resolutions of a numerical model to obtain the best degree of accuracy of numer- 
ical solutions and the best simulation and prediction is an important and urgent 
computational issue. 
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4. Nonlinear Error Growth Dynamics and Predictability 

For a complex nonlinear chaotic system such as the atmosphere or the climate the 
intrinsic randomness in the system sets a theoretical limit to its predictability; see 
Lorenz [90ll92| . Beyond the predictability limit, the system becomes unpredictable. 
In the studies of predictability, the Lyapunov stability theory has been used to 
determine the predictability limit of a nonlinear dynamical system. The Lyapunov 
exponents give a basic measure of the mean divergence or convergence rates of 
nearby trajectories on a strange attractor, and therefore may be used to study the 
mean predictability of chaotic system; see Eckmann and Ruelle [24|, Wolf et al. 
[153] and Fraedrich [28]. Recently, local or finite-time Lyapunov exponents have 
been defined for a prescribed finite-time interval to study the local dynamics on 
an attractor Kazantsev [57], Ziehmann et al [161j and Yoden and Nomura |155j . 
However, the existing local or finite-time Lyapunov exponents, which are same as 
the global Lyapunov exponent, are established on the basis of the fact that the 
initial perturbations are sufficiently small such that the evolution of them can be 
governed approximately by the tangent linear model (TLM) of the nonlinear model, 
which essentially belongs to linear error growth dynamics. Clearly, as long as an 
uncertainty remains infinitesimal in the framework of linear error growth dynamics 
it cannot pose a limit to predictability. To determine the limit of predictability, 
any proposed local Lyapunov exponent must be defined with the respect to the 
nonlinear behaviors of nonlinear dynamical systems Lacarra and Talagrand |67j 
and Mu pT4] . 

In view of the limitations of linear error growth dynamics, it is necessary to 
propose a new approach based on nonlinear error growth dynamics for quantifying 
the predictability of chaotic systems. Ding and Li [22] and Li et. al. [79] have 
presented a nonlinear error growth dynamics which applies fully nonlinear growth 
equations of nonlinear dynamical systems instead of linear approximation to error 
growth equations to discuss the evolution of initial perturbations and employed it 
to study predictability. 

For an ?7-dimensional nonlinear system, the dynamics of small initial perturba- 
tion 6o — 6 (to) £ K" about about an initial point xq = x{to) in the n-dimensional 
phase space are governed by the nonlinear propagator r]{xo, So,t), which propagates 
the initial error forward to the error at the time t = to + t: 

S{to + t) = 77(.to,(5o,t)(5o. 

Then the nonlinear local Lyapunov exponent (NLLE) is defined by 

Ai(xo,do,T) = - In T-^ . 

T \\M\ 

This indicates that Ai(a;o, i5o, t) depends generally on the initial state Xq in the 
phase space, the initial error Sq, and evolution time r. The NLLE is quite dif- 
ferent from the global Lyapunov exponent (GLE) or the local Lyapunov exponent 
(LLE) based on linear error dynamics. If we study the average predictability of the 
whole system, the whole ensemble mean of the NLLE, Xi{5o, t) =< Ai(a;o, 5q, t) > 
should be introduced, where the symbol < • > denotes the ensemble average. Then 
the average predictability limit of a chaotic system could be quantitatively deter- 
mined using the evolution of the mean relative growth of the initial error (RGIE) 
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E{do,T) = exp(Ai ((5o, t)t). According to the chaotic dynamical system theory and 
probabihty theory, the saturation theorem of RGIE [S^ may be obtained as follows. 

Theorem 4.1. [22j For a chaotic dynamic system, the mean relative growth of ini- 
tial error (RGIE) will necessarily reach a saturation value in a finite time interval. 

Once the RGIE reaches the saturation, at the moment almost all predictability of 
chaotic dynamic systems is lost. Therefore, the predictability limit can be defined 
as the time at which the RGIE reaches its saturation level. 

If the first NLLE Ai(a;o, Sq, t) along the most rapidly growing direction has been 
obtained, for an n-dimensional nonlinear dynamic system the first m NLLE spectra 
along other orthogonal directions can be successively determined by the growth 
rate of the volume Vm of an m-dimensional subspacc spanned by the m initial error 
vectors 6m{to) = {Si{to), ■ ■ ■ ,^m(io): 

for m = 2, 3, - • • ,7i, where A,; = Ai(a;o, (Soj ''") is the i-th NLLE of the dynamical 
system. Correspondingly the whole ensemble mean of the i-th NLLE is defined as 
Ai(<5o,T) =< Ai(a;o,(5o,r) >. 

For a chaotic system, each error vector tends to fall along the local direction of 
the most rapid growth. Due to the finite precision of computer, the collapse toward 
a common direction causes the orientation of all error vectors to become indistin- 
guishable. This problem can be overcome by the repeated use of the Gram-Schmidt 
reorthogonalization (GSR) procedure on the vector frame. Giving a set of error vec- 
tors • • • , <5„}, the GSR provides the following orthogonal set {5'i, • • • , 5'^}: 

S'l = <5i, 



The growth rate of the m-dimensional volume can be calculated by the use of the 
first TO orthogonal error vectors, and then the first m NLLE spectra can be obtained 
correspondingly. 

On the other hand, we introduce the local ensemble mean of the NLLE in or- 
der to measure predictability of specified state with certain initial uncertainties in 
the phase space and to investigate distribution of predictability limit in the phase 
space |22j . Assuming that all initial perturbations with the amplitude and random 
directions are on a spherical surface centered at a initial point xq, the local en- 
semble mean of the NLLE relative to Xq is defined as Ai(xo, t) ~< A(xo, e, t) >ff 
for N ^ oo, and then the local average predictability limit of a chaotic system 
at the point xq could be quantitatively determined by examining the evolution of 
the local mean relative growth of initial error (LRGIE) E{xo,t) = exp(AL(a;o, t)t). 
The local ensemble mean of the NLLE different from the whole ensemble mean 
of the NLLE could show local error growth dynamics of subspace on an attrac- 
tor in the phase space. Moreover, in practice the local average predictability limit 
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itself might be regarded as a predictand to provide an estimation of accuracy of 
prediction results. 

The nonlinear error growth theory mentioned above provides a new idea for 
predictability study. However, a great deal of work, including the theory itself, is 
needed. For a real system such as the atmosphere and ocean, the further studies 
related to the following questions are needed: 

(1) Quantitative estimates of the temporal-spatial characteristics of the pre- 
dictability limit of different variables of the atmosphere and ocean by use 
of observation data (Chen et al. |14] made a preliminary attempt to this 
aspect). 

(2) Relationships among the predictability limits of motion on various time and 
space scales. 

(3) Disclosure of the mechanisms influencing predictability from the view of 
nonlinear error growth dynamics. 

(4) Predictability limit varying with changes of initial perturbations. 

(5) Decadal change of the predictability limit. 

(6) Predictability of extreme events. 

(7) Prediction of predictability. 

5. Climate Variability and Successive Bifurcation 

Understanding climate variability and related physical mechanisms and its ap- 
plications to climate prediction and projection are the primary goals in the study of 
climate dynamics. One of problems of climate variability research is to understand 
and predict the periodic, quasi-periodic, aperiodic, and fully turbulent characteris- 
tics of large-scale atmospheric and oceanic flows. Bifurcation theory enables one to 
determine how qualitatively different flow regimes appear and disappear as control 
parameters vary; it provides us, therefore, with an important method to explore 
the theoretical limits of predicting these flow regimes. 

For this purpose, the ideas of dynamical systems theory and nonlinear functional 
analysis have been applied so far to climate dynamics mainly by careful numerical 
studies. These were pioneered by Lorenz |90l [U, Stommel [139j . and Veronis 
[14411146] among others, who explored the bifurcation structure of low-order models 
of atmospheric and oceanic flows. 

Recently, pseudo-arclength continuation methods have been applied to atmo- 
spheric (Legras and Ghil [70]) and oceanic (Speich et al. |137| and Dijkstra pT] ) 
models with increasing horizontal resolution. These numerical bifurcation studies 
have produced so far fairly reliable results for two classes of geophysical flows: (i) 
atmospheric flows in a periodic mid-latitude channel, in the presence of bottom 
topography and a forcing jet; and (ii) oceanic flows in a rectangular mid-latitude 
basin, subject to wind stress on its upper surface; see among others Charney and 
DeVore [E], Pedlosky [H^, Legras and Ghil [70] and Jin and Ghil [53] for saddle- 
node and Hopf bifurcations in the the atmospheric channel, and [llll [^ H^I127l[TTni 
linTllSnKSKnSldlldlS] for saddlc-nodc, pitchfork or Hopf or global bifurcation 
in the oceanic basin. 

Apparently, further numerical bifurcation studies are inevitably necessary. Typi- 
cal problems include 1) continuation algorithms (pseudo-arclength methods and sta- 
bility analysis) applied to large-dimensional dynamical systems (discretized PDEs), 
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2) Galerkin approach using finite-element discretization together with 'homotopic' 
meshes that can deform continuously from a domain to another. 

Another important further direction is to rigorously conduct bifurcation and 
stability analysis for the original partial differential equations models associated 
with typical phenomena. Some progresses have been made in this direction; see 
among others [46l|45l [471 E] . It is clear that much more effort is needed; see also 
Section 6 below. 

Furthermore, very little is known for theoretical and numerical investigations 
on the bifurcations of coupled systems, which are of practical significance for the 
coupled dynamics. It is also practically important to study the variability and 
dynamics of systems under varying external forcing, e.g., the features, processes 
and dynamics of weather and climate varies with the global warming. 

Hereafter in next few subsections, we present some issues on a few specific phys- 
ical phenomena. 

5.1. Wind-driven and thermohaline circulations. For the ocean, basin-scale 
motion is dominated by wind-driven (horizontal) and thermohaline (vertical) cir- 
culations. Their variability, independently and interactively, may play a significant 
role in climate changes, past and future. The wind-driven circulation plays a role 
mostly in the oceans' subannual-to-interannual variability, while the thermohaline 
circulation is most important in decadal-to-millenial variability. 

The thermohaline circulation (THC) is highly nonlinear due to the combined 
effects of the temperature and the salinity on density (Meinckel et. at. |113j : 
Rahmstorf |124| ). which cause the existence of multiple equilibria and thresholds in 
the THC. The abrupt climate is related to the shift between multiple equilibria flow 
regimes in the THC. The sensitivity of the THC to anthropogenic climate forcing 
is still an open question (Rahmstorf [124j : Meincke et. al. [113j ; Thorpe et al., 
[142] ). THis is closely related to the question on whether an abrupt breakdown of 
the THC can result from global warming. In particular, there are two different but 
connected the stability and transitions associated with the problem. The first is 
stability and transitions of the solutions of the partial differential equation models 
in the phase space, and the second is the structure of the solutions and its transi- 
tions in the physical spaces. Issues related to these transitions appear to be very 
important. One such example is the western boundary current separation. The 
physics of the separation of western boundary currents is a long standing problem 
in physical oceanography. The Gulf Stream in the North Atlantic and Kuroshio 
in the North Pacific have a fairly similar behavior with separation from the coast 
occurring at or close to a fixed latitude. The Agulhas Current in the Indian Ocean, 
however, behaves differently by showing a retrofiection accompanied by ring forma- 
tion. The current rushes southward along the east coast of the African, overshoots 
the southern latitude of this continent and then suddenly it turns eastward and 
flows backward into the Indian Ocean. The North Brazil Current in the equatorial 
Atlantic shows a similar, but weaker, retroflection. 

Mathematically speaking, the boundary-layer separation problem is crucial for 
understanding the transition to turbulence and stability properties of fluid flows. 
This problem is also closely linked to structural and dynamical bifurcation of the 
flow through a topological change of its spatial and phase-space structure. This 
program of research has been initiated by Tian Ma and one of the authors in this 
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article, in collaboration in part with Michael Ghil; see [TOTl [Ml [Ml [Ml ESI [Ml ESI 
[SlllSS]- A great deal of further studies in this direction are needed. 

5.2. Intraseasonal Oscillations (ISOs). Another important source of variability 
is related to ISOs such as the Madden- Julian Oscillation (MJO). MJO is a large- 
scale oscillation (wave) in the equatorial region (Madden and Julian [106|, I107j ) 
and is the dominant component of the intraseasonal (30-90 days) variability in the 
tropical atmosphere (Zhang [160] 1. Although some theories and hypotheses have 
been proposed to understand the MJO, a completely satisfactory dynamical theory 
for the MJO has not yet been established (Holton [44]). 

The basic governing equations used to theoretically analyze and simulate the 
MJO could be found in Wang |148| . From the mathematical point of view, the 
well-posedness, asymptotic behavior of the equations, and bifurcation and stability 
analysis arc the first theoretical questions to be answered. 

The observation diagnosis indicates that there is a rich multiscale structure of the 
MJO, and scale interactions might play an important role in the MJO (Slingo et al. 
[136] ). However, whether the scale interactions are essential for the scale selection of 
the MJO is an important open question (Wang [14811160] ). It is therefore necessary 
to develop a multiscale analysis theory of the multiscale model to study the upscale 
energy transfer and to recognize the formation of large-scale structure of the system 
through multiscale interactions. Two basic aspects might be involved into this area. 
One is the significance of short-term cycle in the life cycle of the inherent large-scale 
structure of a dynamical system, e.g., the diurnal cycle to the MJO. The other is 
how the mesoscale and synoptic-scale systems go through certain organized action 
or stochastic dynamics to form a massive behavior and influence movement of this 
large-scale structure. 

Recently, extratropical ISO or mid-high latitude low-frequency variability (LFV) 
has been revealed, e.g., the 70-day oscillation found over the North Atlantic (Felik, 
Gihl and Simonnet, US]; Keppenne, Marcus, Kimoto and Ghil [5H]; Lau, Sheu 
and Kang [69]) that is related to LFV of North Atlantic Oscillation (NAO). The 
dynamics of extratropical ISO or LFV is clearly an attractive field in the future. 

5.3. ENSO. The El Nino-Southern Oscillation (ENSO) is the known strongest 
interannual climate variability associated with strong atmosphere-ocean coupling, 
which has significant impacts on global climate. ENSO is in fact a phenomenon 
that warm events (El Nifio phase) and clod events (La Nina phase) in the equato- 
rial eastern Pacific SST anomaly occur by turns, which associated with persistent 
weakening or strengthening in the trade winds by turns. 

It is convenient, effective, and easily understandable to employ the simplified 
coupled dynamical models to investigate some essential behaviors of ENSO dynam- 
ics. The simplest and leading theoretical models for ENSO are the delayed oscillator 
model (Schopf and Suarez [132| . Battisti and Hirst |4]) and the recharge oscilla- 
tor model (Jin |52]). However, the basic shortcoming of those highly simplified 
models cannot account for the observed irregularity of ENSO, although they could 
qualitatively explain the average features of an ENSO cycle. Hence further study 
is inevitably necessary. Another simplified coupled ocean-atmosphere model which 
can be used to predict ENSO event is the Zcbiak and Cane (ZC) model (Zebiak and 
Cane [T58] ). The atmosphere model is the Gill-type (Gill [41]; Neehn et al. [TT6] ). 
and the ocean model consists of a shallow-water layer with an embedded mixed 
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layer. Also, we would like to mention the intriguing behavior of Boolean Delay 
Equations (BDE) in the ENSO context; see Ghil, Zaliapin and Coluzzi [40j and the 
references therein. It is worth studying the dynamical bifurcation and stability of 
solutions of this kind of simplified models to understand the phase transitions of 
ENSO and its low-frequency variability; see also the dynamical bifurcation theory 
in the next section. 

The ENSO coupling processes and dynamics under the global warming by using 
successive bifurcation theory and the predictability of ENSO by using of nonlinear 
error growth theory are also of considerable practical importance. However, an 
interesting current debate is whether ENSO is best modeled as a stochastic or 
chaotic system - linear and noise-forced, or nonlinear oscillatory and unstable? It 
is obvious that a careful fundamental level examination of the problem is crucial. 

6. New Dynamical Systems Theories and Geophysical Applications 

6.1. Introduction and motivation. As mentioned earlier, most studies on bifur- 
cation issues in geophysical fluid dynamics so far have only considered systems of 
ordinary differential equations (ODEs) that are obtained by projecting the PDEs 
onto a finite-dimensional solution space, either by finite differencing or by truncat- 
ing a Fourier expansion (see Ghil and Childress [32] and further references there). 

A challenge mathematical problem is to conduct rigorous bifurcation and sta- 
bility analysis for the original partial differential equations (PDEs) that govern 
geophysical flows. Progresses in this area should allow us to overcome some of the 
inherent limitations of the numerical bifurcation results that dominate the climate 
dynamics literature up to this point, and to capture the essential dynamics of the 
governing PDE systems. 

Recently, Ma and Wang initiated a study on a new dynamic bifurcation and 
stability theory for dynamical systems. This bifurcation theory is centered at a 
new notion of bifurcation, called attractor bifurcation for dynamical systems, both 
finite dimensional and infinite dimensional. The main ingredients of the theory 
include a) the attractor bifurcation theory, b) steady state bifurcation for a class 
of nonlinear problems with even order non-degenerate nonlinearities, regardless of 
the multiplicity of the eigenvalues, and c) new strategies for the Lyapunov-Schmidt 
reduction and the center manifold reduction procedures. The general philosophy 
is that we first derive general existence of bifurcation to attractors, and then we 
classify the bifurcated attractors to derive detailed dynamics including for instance 
stability of the bifurcated solutions. 

The bifurcation theory has been applied to various problems from science and en- 
gineering, including, in particular, the Kuramoto-Sivashinshy equation, the Cahn- 
Hillard equation, the Ginzburg-Landau equation, Reaction-Diffusion equations in 
Biology and Chemistry, and the Benard convection problem and the Taylor problem 
in classical fluid mechanics; see a recent book by Ma and Wang |100j . For appli- 
cations to geophysical fluid dynamics problems, we have carried out the detailed 
bifurcation and stability analysis for 1) the stratified Boussincsq equations [17], 2) 
the doubly-diffusive modes (both 2D and 3D) [HHS]. 

We proceed with a simple example to illustrate the basic motivation and ideas 
behind the attractor bifurcation theory. For x ~ {xi,X2) € K^, the system x ~ 
Xx— (xf , a;!) -|-o(|x|^) bifurcates from {x, A) = (0, 0) to an attractor T,\ = S^. This 
bifurcated attractor is as shown in Figure 16.11 and contains exactly 4 nodes (the 
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points a, b, c, and d), 4 saddles (the points e, f, g, h), and orbits connecting these 
8 points. From the physical transition point of view, as A crosses 0, the new state 
after the system undergoes a transition is represented by the whole bifurcated 
attractor Ea, rather than any of the steady states or any of the connecting orbits. 
The connecting orbits represents transient states. Note that the global attractor 
is the 2D region enclosed by E^. Wc point out here that the bifurcated attractor 
is different from the study on global attractors of a dissipative dynamical system- 
both finite and infinite dimensional. Global attractor studies the global long time 
dynamics (see among others [571 [U [T71 (TS] ) , while the bifurcated attractor provides 
a natural object for studying dynamical transitions [lOOi 11051 1103j . 
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Figure 6.1. A bifurcated attractor containing 4 nodes (the points 
a, b, c, and d), 4 saddles (the points e, f, g, h), and orbits connect- 
ing these 8 points. 

One important characteristic of the new attractor bifurcation theory is related 
to the asymptotic stability of the bifurcated solutions. This characteristic can be 
viewed as follows. First, as an attractor itself, the bifurcated attractor has a basin of 
attractor, and consequently is a useful object to describe local transitions. Second, 
with detailed classification of the solutions in the bifurcated attractor, we are able 
to access not only the asymptotic stability of the bifurcated attractor, but also the 
stability of different solutions in the bifurcated attractor, providing a more complete 
understanding of the transitions of the physical system as the system parameter 
varies. Third, as Kirchgassner [59j indicated, an ideal stability theorem would 
include all physically meaningful perturbations and today we are still far from this 
goal. In addition, fluid flows are normally time dependent. Therefore bifurcation 
analysis for steady state problems provides in general only partial answers to the 
problem, and is not enough for solving the stability problem. Hence from the 
physical point of view, attractor bifurcation provides a nature tool for studying 
transitions for deterministic systems. 

6.2. A brief account of the attractor bifurcation theory. We now briefly 
present this new attractor bifurcation theory and refer the interested readers to 
|100| I105j for details of the theory and its various applications. 
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We start with a basic state (p of the system, a steady state solution of (|3.1|1 
Then consider Lp = u + Cp. Then (|3.ip becomes 



utl 

(6.1) —^Lxu + G{u,X), 

(6.2) u(0) = uo, 

where i? and i?i are two Hilbert spaces such that Hi H he a, dense and compact 
inclusion, 

The mapping u : [0, oo) — > _ff is the unknown function, A S R is the system 
parameter, and L\ : Hi — > i7 is a family of linear completely continuous fields 
depending continuously on A G M, such that 

Lx = —A + B\ a sectorial operator, 

(6.3) A : Hi ^ H a linear homeomorphism, 
B\ : Hi H a linear compact operator. 

It is known that L\ generates an analytic semigroup {e~*^^}t>o and we can define 
fractional power operators L" for a e M with domain Ha = D{L") such that 
Hai C Ha2 is compact if ai > a2, H ~ Hq and Hi = Ha=i- 

Furthermore, we assume that for some 6 < 1 the nonlinear operator G(-,A) : 
He Hq is a C bounded operator (?- > 1), and 

(6.4) G(u,A) = odiulle), V A e K. 

Definition 6.1 (Ma & Wang [98l [100]). (1) We say that the equation fO]) bi- 
furcates from (u. A) = (0, Ao) to an invariant set Sa, if there exists a se- 
quence of invariant sets of i6.1\) . such that ^ and 

lim A„ = Ao, lim max ||a;|| = 0. 

(2) // the invariant sets Sa are attractors of h6.1\) . then the bifurcation is called 
an attractor bifurcation. 




Figure 6.2. Definition of attractor bifurcation. 
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Let {/3fc(A) e C I fc = 1, 2, • • • } be the eigenvalues of Lx = —A + Bx (counting 
the multiphcities). Suppose that 







f <0 


if A < Ao 


(6.5) 


i?eA(A) 1 


' =0 


if A = Ao 






i>„ 


if A > Ao 



(6.6) Re(3j{Xo) < for j>m+ 1. 

Let-Eo be the eigenspace of Lx at Ao 

Eq^ [j {ueHl (Lao-A(Ao))^u = o, fceN}. 

l<2<m 

By (|6.5p we know that dimi^o = "i. 

In physical terms, the above properties of eigenvalues are called principle of 
exchange of stabilities (PES). They can be verified in many physical systems. It 
is obvious that they are necessary for linear instability and the following attractor 
bifurcation theorem demonstrates that they lead to bifurcation to an attractor. 

Theorem 6.2 (Ma & Wang [100l[98])- Assume that ifO) ) and ifO)) hold true, and 
let u = be a locally asymptotically stable equilibrium point of 16.1]) at X = Xq. 
Then we have the following assertions. 

(1) Equation id. 1\) bifurcates from {u, X) = (0, Aq) to an attractor T,x for A > Aq 
with m — 1 < dim^A < m, which is an (m — \)- dimensional homological 
sphere. 

(2) For any ux G Sa; u\ can be expressed as 

ux = vx + o(||t;A||), 

where vx (z Eq. 

(3) There exists a neighborhood U d H of u ^ such that Ea attracts U \ T, 
where T is the stable manifold of u = with codimension m. In particular, 
if u ~ is global asymptotically stable for i6.1]) at X ~ Xq and 116.1]) has a 
global attractor for any X near Xq, then Sa attracts H \ T. 

With this general theorem, along with other important ingredients addressed 
in Ma and Wang [1001 1105] . in our disposal, the main strategy to conduct the 
bifurcation analysis consists of 1) the existence of attractor bifurcation, and 2) 
classification of solutions in the bifurcated attractor. We refer the interested readers 
to the books |100|, I105j for details; see also the discussion below on the Raylcigh- 
Bcnard convection. 

6.3. Dynamic bifurcation in classical fluid dynamics. The new dynamic bi- 
furcation theory has been naturally applied to problems in fluid dynamics, including 
the Rayleigh-Benard convection, the Taylor problem, and the parallel shear flows, 
by Ma and Wang and their collaborators. The study of these basic problems on 
the one hand plays an important role in understanding the turbulent behavior of 
fluid flows, and on the other hand often leads to new insights and methods toward 
solutions of other problems in sciences and engineering. 

To illustrate the main ideas of the applications of the theory, we consider now the 
Rayleigh-Benard convection problem. Linear theory of the Raylcigh-Benard prob- 
lem were essentially derived by physicists; see, among others, Chandrasekhar [TO] 
and Drazin and Reid [23] . Bifurcating solutions of the nonlinear problem were first 
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constructed formally by Malkus and Veronis [109j . The first rigorous proofs of the 
existence of bifurcating solutions were given by Yudovich [156| I157j and Rabinowitz 
|123| . Yudovich proved the existence of bifurcating solutions by a topological de- 
gree argument. Earlier, however, Velte [143] had proved the existence of branching 
solutions of the Taylor problem by a topological degree argument as well. The have 
been many studies since the above mentioned early works. However, as we shall 
see, a rigorous and complete understanding of the problem is not available until the 
recent work by Ma and Wang using the attractor bifurcation theory. 

To illustrate the ideas, we start by recalling the basic set-up of the problem. Let 
the Rayleigh number be 

^ _ ga{% - fi)h^ 

The Benard convection is modeled by the Boussinesq equations. In their nondi- 
mensional form, these equations are written as follows: 

1 



Pr 

(6.7) dT 

'dt 

div u = 



du , „ 
— + (u • V)m + Vp 



- Alt - y/RTk = 0, 
+ {u ■ V)T - y/R.U3 - AT = 0, 



Here the unknown functions are {u,T,p), which are the deviations from the basic 
field. The non-dimensional domain is = Z) x (0, 1) C M'^, where D C is 
an open set. The coordinate system is given hy x = {xi,X2,X3) G M.^. They are 
supplemented with the following initial value conditions 

(6.8) (u,T) = (uo,To) &tt = 0. 

Boundary conditions are needed at the top and bottom and at the lateral bound- 
ary dD X (0, 1). At the top and bottom boundary {x^ ~ 0, 1), either the so-called 
rigid or free boundary conditions are given 

T = 0, M = (rigid boundary), 

T = 0, W3 = 0, — = (free boundary). 

0x3 

Different combinations of top and bottom boundary conditions are normally used 
in different physical setting such as rigid-rigid, rigid-free, free-rigid, and free-free. 

On the lateral boundary dD x [0, 1], one of the following boundary conditions 
are usually used: 

(1) Periodic condition: 

(6.10) {u,T){xi + kiLi,X2 + k2L2,X3) = {u,T){xi,X2,X3), 

for any ki,k2 G Z. 

(2) Dirichlet boundary condition: 

dT 

(6.11) u = 0, r = (or — = 0); 

on 

(3) Free boundary condition: 

(6.12) r = o, u„ = 0, ^=0, 
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where n and r are the unit normal and tangent vectors on dD x [0, 1] 
respectively, and Un = u ■ n, Ut = u ■ t . 
By using the attractor bifurcation theory, the following results have been ob- 
tained by Ma and Wang in [93 [102 [lOOl [M] • 

(1) When the Rayleigh number R crosses the first critical Rayleigli number 
the Rayleigh-Benard problem bifurcates from the basic state to an attractor 
An, homologic to S'™~^, where m is the multiplicity of Rc as an eigenvalue 
of the linearized problem near the basic solution, for all physically sound 
boundary conditions, regardless of the geometry of the domain and the 
multiplicity of the eigenvalue Rc for the linear problem 

(2) Consider the 3D Benard convection in = (0, Li) x (0, L2) x (0,1) with 
free top-bottom and periodic horizontal boundary conditions, and with 

P P 1 

(6.13) + 7T ~ o some ki,k2 <E Z. 
Then 

(6.14) Ar . 



if L2 = Vfc2 - \Li, fc = 2,3, •■ 
otherwise. 



(3) For the 3D Benard convection in Q. ~ (0, L)^ x (0, 1) with free boundary 
conditions and with 

2 - 2I/3 



(6.15) 0<^'<2i/3_i- 

The bifurcated attractor Ar consists of exactly eight singular points and 
eight heteroclinic orbits connecting the singular points, as shown in Figure 
16.11 with 4 of them being minimal attractors, and the other 4 saddle points. 
The proof of these results is carried out in the following steps: 

1) It is classical to put the Boussinesq equations (|6.7p in the abstract form 
(|6.ip . Then it is easy to see that the linearized operator is a self-adjoint 
operator., and both conditions (|6.3p and (|6.4p arc satisfied. 

2) As the linearized operator is symmetric, it is not hard to verify the PES 
(|6.5p and (|6.6p for the Boussinesq equations. 

3) To derive a general attractor theorem for the Benard convection regardless 
of the domain and the boundary conditions, we need to verify the asymp- 
totic stability of the basic solution [u, T) — Q &t the critical Reynolds 
number. This can be achieved by a general stability theorem, derived by 
Ma and Wang in [HZ] . 

4) The detailed structure of the solutions in the bifurcated attractor can be 
derived using a new approximation formula for the center manifold function; 
see Ma and Wang [lOOl [105] . 

We remark here that the high dimensional sphere S*"'"^ contains not only the 
steady states generated by symmetry groups inherited in the problem, but also 
many transient states, which were completely missed by any classical theories. 
These transient states are highly relevant in geophysical fluid dynamics and climate 
dynamics. We believe the climate low frequency climate variabilities discussed in 
the previous section are related to certain transient states, and further exploration 
in this direction are certainly important and necessary. 
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6.4. Dynamic bifurcation and stability in geophysical fluid dynamics. As 

mentioned earlier, the theory has been applied to models in geophysical fluid dy- 
namics models including the doubly-diffusive models, and rotating Boussinesq equa- 
tions. We now briefly address these applications in turn. 

Rotating Boussinesq Equations: Rotating Boussinesq equations are basic 
model in atmosphere and ocean dynamical models. In |47] . for the case where the 
Prandtl number is greater than one, a complete stability and bifurcation analysis 
near the first critical Rayleigh number is carried out. Second, for the case where 
the Prandtl number is smaller than one, the onset of the Hopf bifurcation near the 
first critical Rayleigh number is established, leading to the existence of nontrivial 
periodic solutions. The analysis is based on a newly developed bifurcation and 
stability theory for nonlinear dynamical systems as mentioned above. 

Double-Diffusive Ocean Model: Double-diffusion was first originally dis- 
covered in the 1857 by Jevons [51j . forgotten, and then rediscovered as an "occano- 
graphic curiosity" a century later; see among others Stommel, Arons and Blanchard 
|140j . Veronis [145j . and Baines and Gill [3]. In addition to its effects on oceanic 
circulation, double-diffusion convection has wide applications to such diverse fields 
as growing crystals, the dynamics of magma chambers and convection in the sun. 

The best known doubly-diffusive instabilities are "salt-fingers" as discussed in 
the pioneering work by Stern [138| . These arise when hot salty water lies over cold 
fresh water of a higher density and consist of long fingers of rising and sinking water. 
A blob of hot salty water which finds itself surrounded by cold fresh water rapidly 
loses its heat while retaining its salt due to the very different rates of diffusion 
of heat and salt. The blob becomes cold and salty and hence denser than the 
surrounding fluid. This tends to make the blob sink further, drawing down more 
hot salty water from above giving rise to sinking fingers of fiuid. 

In [46( I45j , we present a bifurcation and stability analysis on the doubly-diffusive 
convection. The main objective is to study 1) the mechanism of the saddle- node 
bifurcation and hysteresis for the problem, 2) the formation, stability and transi- 
tions of the typical convection structures, and 3) the stability of solutions. It is 
proved in particular that there are two different types of transitions: continuous 
and jump, which are determined explicitly using some physical relevant nondimen- 
sional parameters. It is also proved that the jump transition always leads to the 
existence of a saddle-node bifurcation and hysteresis phenomena. 

However, there are many issues still open, including in particular to use the 
dynamical systems tools developed to study some of the issues raised in the previous 
sections. 

6.5. Stability and transitions of geophysical flows in the physical space. 

Another important area of studies in geophysical fluid dynamics is to study the 
structure and its stability and transitions of flows in the physical spaces. 

A method to study these important problems in geophysical fluid dynamics is 
a recently developed geometric theory for incompressible flows by Ma and Wang 
[101] . This theory consists of research in directions: 1) the study of the structure 
and its transitions/evolutions of divergence-free vector fields, and 2) the study of 
the structure and its transitions of velocity fields for 2-D incompressible fluid flows 
governed by the Navier-Stokes equations or the Euler equations. The study in 
the first direction is more kinematic in nature, and the results and methods devel- 
oped can naturally be applied to other problems of mathematical physics involving 
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divergence-free vector fields. In fluid dynamics context, the study in the second 
direction involves specific connections between the solutions of the Navier-Stokes 
or the Euler equations and flow structure in the physical space. In other words, 
this area of research links the kinematics to the dynamics of fluid flows. This is 
unquestionably an important and difficult problem. 

Progresses have been made in several directions. First, a new rigorous char- 
acterization of boundary layer separations for 2-D viscous incompressible flows is 
developed recently by Ma and Wang, in collaboration in part with Michael Ghil 
[101] . The nature of flow's boundary layer separation from the boundary plays a 
fundamental role in many physical problems, and often determines the nature of 
the flow in the interior as well. The main objective of this section is to present a 
rigorous characterization of the boundary layer separations of 2D incompressible 
fluid flows. This is a long standing problem in fluid mechanics going back to the pi- 
oneering work of Prandtl [121j in 1904. No known theorem, which can be applied to 
determine the separation, is available until the recent work by Ghil, Ma and Wang 
[36l[37j, and Ma and Wang [93l[96], which provides a first rigorous characterization. 
Interior separations arc studied rigorously by Ma and Wang [99]. The results for 
both the interior and boundary layer separations arc used in Ma and Wang [102j for 
the transitions of the Couctte-Poiseuille and the Taylor-Couette-Poiseuille flows. 

Another example in this area is the justification of the roll structure (e.g. rolls) in 
the physical space in the Rayleigh-Benard convection by Ma and Wang [97[|104| . We 
note that a special structure with rolls separated by a cross channel flow derived in 
[104j has not been rigorously examined in the Benard convection setting although it 
has been observed in other physical contexts such as the Branstator-Kushnir waves 
in the atmospheric dynamics [T] [S^ . 

With this theory in our disposal, the structure/patterns and their stability and 
transitions in the underlying physical spaces for those problems in fluid dynamics 
and in geophysical fluid dynamics can be classified. In particular, this theory has 
been used to study the formation, persistence and transitions of flow structures 
including boundary layer separation including the Gulf separation, the Hadley cir- 
culation, and the Walker circulation. Further investigation appears to be utterly 
important and necessary. In addition, it appears that such theoretical and numer- 
ical studies will lead to better predications on weather and climate regimes. 
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